/******************************************************************************
Safer in School? The Impact of Compulsory Schooling on Maltreatment and Associated Harms
by Adam A. Dzulkipli, Nicole Black, David W. Johnston, and Leonie Segal

Description: This program recreates all of the main figures (i.e. Figures 1 to 3)
presented in the main text of the paper. The figures are exported as individual
PDF files, and can then be combined together manually using a typesetting
software such as LaTeX or Microsoft Word.

Note that directory paths have been removed to preserve privacy. Please insert 
the relevant directories in the global macros under "Setup" before running the
.do file. This code was written and implemented in Stata MP 19.5, with
a computational time of approximately <1 minute. 

Packages required include:
•	ftools
•	reghdfe
•	estout
•	rdrobust
•	pzms
•	rddisttestk (.ado file by Brigham Frandsen available from
	https://sites.google.com/view/brighamfrandsen/software?authuser=0)
•	rddensity
•	rdhonest

*****************************************************************************/

********************************************************************************
*	Setup
********************************************************************************
clear all
set more off
macro drop all
capture log close

graph set window fontface "Times New Roman"
set scheme s2color

global raw
global tmp
global clean
global output

********************************************************************************
*	Figure 1: RD plots for school enrolment
********************************************************************************

use "${clean}analysis_final.dta", clear

* Restrict sample to children for whom we have schooling records
keep if rec_census15 == 1 

/*
  Calculate school enrolment rates by day-of-birth relative the reform cut-off
  (i.e. 1 January 1993), grouped together in monthly bins
*/
local averages enrol16 enrol17

foreach var in `averages'{
	bysort h1: egen `var'_mean = mean(`var')
	
	sort h1 h2
	bysort h1 (h2): replace `var'_mean = . if  _n != 1
}

foreach var in `averages'{
	bysort cps_age16_binary h1: egen `var'_mean_cp = mean(`var')
	
	sort cps_age16_binary h1 h2
	bysort cps_age16_binary h1 (h2): replace `var'_mean_cp = . if  _n != 1
}

* Define macros and graph settings
global bw1 "h2 >= -721 & h2 <= 721"
global bwr1 "h2 >= 0 & h2 <= 721"
global bwl1 "h2 < 0 & h2 >= -721"
global xtitle "Days born before/after January 1993"

global cp0 "cps_age16_binary == 0"
global cp1 "cps_age16_binary == 1"
global xlabel "xlab(-720(240)720, nogrid labsize(4.5)) xscale(range(-780 780))"
global graph_opt "xline(0, lwidth(thin) lpattern(dash) lcolor(black)) legend(off) graphregion(color(white))"

*** Panel A: RD plots for school enrolment (entire sample)
* Enrolment at age 16
global yvar "enrol16"
global ylabel "ylab(0.7(0.1)1, nogrid labsize(4.5)) ytitle("")"
global title "Enrolled in School at Age 16"

twoway (scatter ${yvar}_mean h2 if $bw1, msymbol(triangle) mcolor(gs10) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1, lcolor(black)), ///
	title($title, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt

graph export "${output}Figure1_a.pdf", replace

* Enrolment at age 17
global yvar "enrol17"
global ylabel "ylab(0.5(0.1)1, nogrid labsize(4.5))"
global title "Enrolled in School at Age 17"

twoway (scatter ${yvar}_mean h2 if $bw1, msymbol(triangle) mcolor(gs10) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1, lcolor(black)), ///
	title($title, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt
	
graph export "${output}Figure1_b.pdf", replace

*** Panel B: RD plots for school enrolment (no past CPS involvement)
* Enrolment at age 16
global yvar "enrol16"
global ylabel "ylab(0.7(0.1)1, nogrid labsize(4.5)) ytitle("")"
global title0 "Enrolled in School at Age 16 (Non-CPS)"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp0, msymbol(diamond) mcolor(gs7) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp0, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp0, lcolor(black)), ///
	title($title0, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt

graph export "${output}Figure1_c.pdf", replace

* Enrolment at age 17
global yvar "enrol17"
global ylabel "ylab(0.5(0.1)1, nogrid labsize(4.5))"
global title0 "Enrolled in School at Age 17 (Non-CPS)"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp0, msymbol(diamond) mcolor(gs7) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp0, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp0, lcolor(black)), ///
	title($title0, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt
	
graph export "${output}Figure1_d.pdf", replace

*** Panel C: RD plots for school enrolment (past CPS involvement)
* Enrolment at age 16
global yvar "enrol16"
global ylabel "ylab(0.7(0.1)1, nogrid labsize(4.5)) ytitle("")"
global title1 "Enrolled in School at Age 16 (Past CPS)"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp1, mcolor(gs4) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp1, lcolor(black)), ///
	title($title1, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt

graph export "${output}Figure1_e.pdf", replace

* Enrolment at age 17
global yvar "enrol17"
global ylabel "ylab(0.5(0.1)1, nogrid labsize(4.5))"
global title1 "Enrolled in School at Age 17 (Past CPS)"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp1, mcolor(gs4) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp1, lcolor(black)), ///
	title($title1, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt
	
graph export "${output}Figure1_f.pdf", replace

********************************************************************************
*	Figure 2: RD plots for new CPS involvement
********************************************************************************

use "${clean}analysis_final.dta", clear

/*
  Calculate new CPS notification rates by day-of-birth relative to the reform 
  cut-off (i.e. 1 January 1993), grouped together in monthly bins
*/
local averages cps_after16_binary

foreach var in `averages'{
	bysort cps_age16_binary h1: egen `var'_mean_cp = mean(`var')
	
	sort cps_age16_binary h1 h2
	bysort cps_age16_binary h1 (h2): replace `var'_mean_cp = . if _n != 1
}

foreach var in `averages'{
	bysort h1: egen `var'_mean = mean(`var')
	
	sort h1 h2
	bysort h1 (h2): replace `var'_mean = . if _n != 1
}

* Define macros and graph settings
global bw1 "h2 >= -721 & h2 <= 721"
global bwr1 "h2 >= 0 & h2 <= 721"
global bwl1 "h2 < 0 & h2 >= -721"
global xtitle "Days born before/after January 1993"

global cp0 "cps_age16_binary == 0"
global cp1 "cps_age16_binary == 1"
global xlabel "xlab(-720(240)720, labsize(4.5)) xscale(range(-780 780))"
global graph_opt "xline(0, lwidth(thin) lpattern(dash) lcolor(black)) legend(off) graphregion(color(white))"

*** Panel A: RD plot for new CPS notifications (entire sample)
global yvar "cps_after16_binary"
global ylabel "ylab(0(0.02)0.1, nogrid labsize(4.5))"
global title "Any CPS Notification from Age 16"

twoway (scatter ${yvar}_mean h2 if $bw1, msymbol(triangle) mcolor(gs10) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1, lcolor(black)), ///
	title($title, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt
	
graph export "${output}Figure2_a.pdf", replace

*** Panel B: RD plots for new CPS notifications (by past CPS involvement)
* Plot for children with no past CPS involvement
global yvar "cps_after16_binary"
global ylabel0 "ylab(0(0.01)0.04, nogrid labsize(4.5))"
global title0 "Any CPS Notification from Age 16 (No Past CPS)"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp0, msymbol(diamond) mcolor(gs7) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp0, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp0, lcolor(black)), ///
	title($title0, size(5.5) color(black)) $ylabel0 ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt
	
graph export "${output}Figure2_b.pdf", replace

* Plot for children with past CPS involvement
global yvar "cps_after16_binary"
global ylabel1 "ylab(0(0.1)0.4, nogrid labsize(4.5))"
global title1 "Any CPS Notification from Age 16 (Past CPS)"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp1, mcolor(gs4) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp1, lcolor(black)), ///
	title($title1, size(5.5) color(black)) $ylabel1 ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt
	
graph export "${output}Figure2_c.pdf", replace

********************************************************************************
*	Figure 3: RD plots for ED utilisation
********************************************************************************

use "${clean}analysis_final.dta", clear

/*
  Calculate ED utilisation rates by day-of-birth relative to the reform 
  cut-off (i.e. 1 January 1993), grouped together in monthly bins
*/
local averages evered_17to21 evered16

* Local averages for entire sample
foreach var in `averages'{
	bysort h1: egen `var'_mean = mean(`var')
	
	sort h1 h2
	bysort h1 (h2): replace `var'_mean = . if _n != 1
}

* Local averages by past CPS involvement
foreach var in `averages'{
	bysort cps_age16_binary h1: egen `var'_mean_cp = mean(`var')
	
	sort cps_age16_binary h1 h2
	bysort cps_age16_binary h1 (h2): replace `var'_mean_cp = . if _n != 1
}

* Define macros and graph settings
global bw1 "h2 >= -721 & h2 <= 721"
global bwr1 "h2 >= 0 & h2 <= 721"
global bwl1 "h2 < 0 & h2 >= -721"
global xtitle "Days born before/after January 1993"

global cp0 "cps_age16_binary == 0"
global cp1 "cps_age16_binary == 1"
global xlabel "xlab(-720(240)720, nogrid labsize(4.5) nogrid) xscale(range(-780 780))"
global graph_opt "xline(0, lwidth(thin) lpattern(dash) lcolor(black)) graphregion(color(white))"

*** Panel A: RD plots for entire sample
* Any ED visit at age 16
global yvar "evered16"
global ylabel "ylab(0(0.1)0.4, labsize(4.5) nogrid) ytitle("")"
global title "Any ED Visit at Age 16"

twoway (scatter ${yvar}_mean h2 if $bw1, msymbol(triangle) mcolor(gs10) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1, lcolor(black)), ///
	title($title, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt legend(off)
	
graph export "${output}Figure3_a.pdf", replace
	
* Any ED visit between ages 17 and 21
global yvar "evered_17to21"
global ylabel "ylab(0.2(0.1)0.8, labsize(4.5) nogrid)"
global title "Any ED Visit between Ages 17 and 21"

twoway (scatter ${yvar}_mean h2 if $bw1, msymbol(triangle) mcolor(gs10) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1, lcolor(black)), ///
	title($title, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt legend(off) 
	
graph export "${output}Figure3_b.pdf", replace

*** Panel B: RD plots by past CPS involvement
* Any ED visit at age 16
global yvar "evered16"
global ylabel "ylab(0(0.1)0.4, labsize(4.5) nogrid) ytitle("")"
global title "Any ED Visit at Age 16"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp0, msymbol(diamond) mcolor(gs7) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp0, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp0, lcolor(black))	///
	(scatter ${yvar}_mean_cp h2 if $bw1 & $cp1, mcolor(gs4) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp1, lcolor(black)), ///
	title($title, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt legend(order(1 "Non-CPS" 8 "Past CPS") position(0) bplacement(neast) size(*1))
	
graph export "${output}Figure3_c.pdf", replace
	
* Any ED visit between ages 17 and 21
global yvar "evered_17to21"
global ylabel "ylab(0.2(0.1)0.8, labsize(4.5) nogrid)"
global title "Any ED Visit between Ages 17 and 21"

twoway (scatter ${yvar}_mean_cp h2 if $bw1 & $cp0, msymbol(diamond) mcolor(gs7) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp0, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp0, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp0, lcolor(black))	///
	(scatter ${yvar}_mean_cp h2 if $bw1 & $cp1, mcolor(gs4) msize(medlarge)) ///
	(lfitci ${yvar} h2 if $bwl1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfitci ${yvar} h2 if $bwr1 & $cp1, ciplot(rarea) fcolor(%20) lwidth(vthin)) ///
	(lfit ${yvar} h2 if $bwl1 & $cp1, lcolor(black)) ///
	(lfit ${yvar} h2 if $bwr1 & $cp1, lcolor(black)), ///
	title($title, size(5.5) color(black)) $ylabel ///
	$xlabel xtitle($xtitle, size(5) color(black)) ///
	$graph_opt legend(order(1 "Non-CPS" 8 "Past CPS") position(0) bplacement(neast) size(*1))
	
graph export "${output}Figure3_d.pdf", replace
